#ifndef ARRAY3_H
#define ARRAY3_H

#include <algorithm>
#include "LosTopos/common/array1.h"
#include <cassert>
#include <vector>

namespace LosTopos {

template<class T, class ArrayT=std::vector<T> >
struct Array3
{
    // STL-friendly typedefs
    
    typedef typename ArrayT::iterator iterator;
    typedef typename ArrayT::const_iterator const_iterator;
    typedef typename ArrayT::size_type size_type;
    typedef long difference_type;
    typedef T& reference;
    typedef const T& const_reference;
    typedef T value_type;
    typedef T* pointer;
    typedef const T* const_pointer;
    typedef typename ArrayT::reverse_iterator reverse_iterator;
    typedef typename ArrayT::const_reverse_iterator const_reverse_iterator;
    
    // the actual representation
    
    int ni, nj, nk;
    ArrayT a;
    
    // the interface
    
    Array3(void)
    : ni(0), nj(0), nk(0), a(0)
    {}
    
    Array3(int ni_, int nj_, int nk_)
    : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_)
    { assert(ni_>=0 && nj_>=0 && nk_>=0); }
    
    Array3(int ni_, int nj_, int nk_, ArrayT& a_)
    : ni(ni_), nj(nj_), nk(nk_), a(a_)
    { assert(ni_>=0 && nj_>=0 && nk_>=0); }
    
    Array3(int ni_, int nj_, int nk_, const T& value)
    : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, value)
    { assert(ni_>=0 && nj_>=0 && nk_>=0); }
    
    Array3(int ni_, int nj_, int nk_, const T& value, size_type max_n_)
    : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, value, max_n_)
    { assert(ni_>=0 && nj_>=0 && nk_>=0); }
    
    Array3(int ni_, int nj_, int nk_, T* data_)
    : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, data_)
    { assert(ni_>=0 && nj_>=0 && nk_>=0); }
    
    Array3(int ni_, int nj_, int nk_, T* data_, size_type max_n_)
    : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, data_, max_n_)
    { assert(ni_>=0 && nj_>=0 && nk_>=0); }
    
    ~Array3(void)
    {
#ifndef NDEBUG
        ni=nj=0;
#endif
    }
    
       const T& operator()(int i, int j, int k) const
       {
          assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
          return a[i+ni*(j+nj*k)];
       }
    
       T& operator()(int i, int j, int k)
       {
          assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
          return a[i+ni*(j+nj*k)];
       }
    
  /*  const T& operator()(size_t i, size_t j, size_t k) const
    {
        assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
        return a[i+ni*(j+nj*k)];
    }
    
    T& operator()(size_t i, size_t j, size_t k)
    {
        assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
        return a[i+ni*(j+nj*k)];
    }*/
    
    bool operator==(const Array3<T>& x) const
    { return ni==x.ni && nj==x.nj && nk==x.nk && a==x.a; }
    
    bool operator!=(const Array3<T>& x) const
    { return ni!=x.ni || nj!=x.nj || nk!=x.nk || a!=x.a; }
    
    bool operator<(const Array3<T>& x) const
    {
        if(ni<x.ni) return true; else if(ni>x.ni) return false;
        if(nj<x.nj) return true; else if(nj>x.nj) return false;
        if(nk<x.nk) return true; else if(nk>x.nk) return false;
        return a<x.a;
    }
    
    bool operator>(const Array3<T>& x) const
    {
        if(ni>x.ni) return true; else if(ni<x.ni) return false;
        if(nj>x.nj) return true; else if(nj<x.nj) return false;
        if(nk>x.nk) return true; else if(nk<x.nk) return false;
        return a>x.a;
    }
    
    bool operator<=(const Array3<T>& x) const
    {
        if(ni<x.ni) return true; else if(ni>x.ni) return false;
        if(nj<x.nj) return true; else if(nj>x.nj) return false;
        if(nk<x.nk) return true; else if(nk>x.nk) return false;
        return a<=x.a;
    }
    
    bool operator>=(const Array3<T>& x) const
    {
        if(ni>x.ni) return true; else if(ni<x.ni) return false;
        if(nj>x.nj) return true; else if(nj<x.nj) return false;
        if(nk>x.nk) return true; else if(nk<x.nk) return false;
        return a>=x.a;
    }
    
    void assign(const T& value)
    { a.assign(value); }
    
    void assign(int ni_, int nj_, int nk_, const T& value)
    {
        a.assign(ni_*nj_*nk_, value);
        ni=ni_;
        nj=nj_;
        nk=nk_;
    }
    
    void assign(int ni_, int nj_, int nk_, const T* copydata)
    {
        a.assign(ni_*nj_*nk_, copydata);
        ni=ni_;
        nj=nj_;
        nk=nk_;
    }
    
    const T& at(int i, int j, int k) const
    {
        assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
        return a[i+ni*(j+nj*k)];
    }
    
    T& at(int i, int j, int k)
    {
        assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
        return a[i+ni*(j+nj*k)];
    }
    
    const T& back(void) const
    { 
        assert(a.size());
        return a.back();
    }
    
    T& back(void)
    {
        assert(a.size());
        return a.back();
    }
    
    const_iterator begin(void) const
    { return a.begin(); }
    
    iterator begin(void)
    { return a.begin(); }
    
    size_type capacity(void) const
    { return a.capacity(); }
    
    void clear(void)
    {
        a.clear();
        ni=nj=nk=0;
    }
    
    bool empty(void) const
    { return a.empty(); }
    
    const_iterator end(void) const
    { return a.end(); }
    
    iterator end(void)
    { return a.end(); }
    
    void fill(int ni_, int nj_, int nk_, const T& value)
    {
        a.fill(ni_*nj_*nk_, value);
        ni=ni_;
        nj=nj_;
        nk=nk_;
    }
    
    const T& front(void) const
    {
        assert(a.size());
        return a.front();
    }
    
    T& front(void)
    {
        assert(a.size());
        return a.front();
    }
    
    size_type max_size(void) const
    { return a.max_size(); }
    
    reverse_iterator rbegin(void)
    { return reverse_iterator(end()); }
    
    const_reverse_iterator rbegin(void) const
    { return const_reverse_iterator(end()); }
    
    reverse_iterator rend(void)
    { return reverse_iterator(begin()); }
    
    const_reverse_iterator rend(void) const
    { return const_reverse_iterator(begin()); }
    
    void reserve(int reserve_ni, int reserve_nj, int reserve_nk)
    { a.reserve(reserve_ni*reserve_nj*reserve_nk); }
    
    void resize(int ni_, int nj_, int nk_)
    {
        assert(ni_>=0 && nj_>=0 && nk_>=0);
        a.resize(ni_*nj_*nk_);
        ni=ni_;
        nj=nj_;
        nk=nk_;
    }
    
    void resize(int ni_, int nj_, int nk_, const T& value)
    {
        assert(ni_>=0 && nj_>=0 && nk_>=0);
        a.resize(ni_*nj_*nk_, value);
        ni=ni_;
        nj=nj_;
        nk=nk_;
    }
    
    void set_zero(void)
    { a.set_zero(); }
    
    size_type size(void) const
    { return a.size(); }
    
    void swap(Array3<T>& x)
    {
        std::swap(ni, x.ni);
        std::swap(nj, x.nj);
        std::swap(nk, x.nk);
        a.swap(x.a);
    }
    
    void trim(void)
    { a.trim(); }
};

// some common arrays

typedef Array3<double, Array1<double> > Array3d;
typedef Array3<float, Array1<float> > Array3f;
//typedef Array3<long long, Array1<long long> > Array3ll;
//typedef Array3<unsigned long long, Array1<unsigned long long> > Array3ull;
typedef Array3<int, Array1<int> > Array3i;
typedef Array3<unsigned int, Array1<unsigned int> > Array3ui;
typedef Array3<short, Array1<short> > Array3s;
typedef Array3<unsigned short, Array1<unsigned short> > Array3us;
typedef Array3<char, Array1<char> > Array3c;
typedef Array3<unsigned char, Array1<unsigned char> > Array3uc;

}

#endif
